Όλα τα απαιτούμενα δεδομένα και αρχεία θα τα βρείτε στον ιστότοπο του μαθήματος στο GitHub

1 Πρώτα βήματα

1.1 Εγκατάσταση βιβλιοθήκης

Εγκατάσταση του πρώτου και κύριου πακέτου (και ενός βοηθητικού)
(Αυτό γίνεται μια φορά σε κάθε Η/Υ - εάν έχει εγκατασταθεί δεν χρειάζεται να δοθεί ξανά η εντολή αυτή)

1.2 Φόρτωση βιβλιοθήκης

Φόρτωση του εν λόγω πακέτου

1.3 Εγκατάσταση βιβλιοθηκών

Με την εντολή packages('insert_package_name_here') γίνεται η εγκατάσταση των πακέτων που μας ενδιαφέρουν

Δοκιμάστε το

1.4 Φόρτωση βιβλιοθηκών

Εντολή για φόρτωση των πακέτων

1.5 Working Directory

Με την ακόλουθη εντολή, βλέπουμε σε ποιόν φάκελο δουλεύουμε:

1.6 Αλλαγή του working directory

Ενώ με την εντολή αυτή, αλλάζουμε τον φάκελο στον οποίο δουλεύουμε και αποθηκεύουμε δεδομένα
[π.χ., setwd("E:/") εάν θέλουμε να δουλεύουμε στον σκληρό δίσκο Ε]

1.7 Εισαγωγή αρχείου

Τώρα ας εισάγουμε το αρχείο το οποίο περιέχει τα δεδομένα τα οποία μας ενδιαφέρουν να αναλύσουμε

Όπως αντιλαμβάνεστε, σήμερα θα δουλέψουμε με δεδομένα τα οποία ήδη υπάρχουν εντός της R. Όμως για την εργασία που θα παραδώσετε εσείς, θα χρειαστεί να φορτώσετε δεδομένα από ένα αρχείο excel, το οποίο θα σας δοθεί στο τέλος της άσκησης. Για να το κάνετε αυτό, θα χρειαστεί να χρησιμοποιήσετε την βιβλιοθηκη readxl και την εντολη read_excel1

Ας δούμε τα δεδομένα μας:

## [1] "Sepal.Length" "Sepal.Width"  "Petal.Length" "Petal.Width"  "Species"
Sepal.Length Sepal.Width Petal.Length Petal.Width Species
5.1 3.5 1.4 0.2 setosa
4.9 3.0 1.4 0.2 setosa
4.7 3.2 1.3 0.2 setosa
4.6 3.1 1.5 0.2 setosa
5.0 3.6 1.4 0.2 setosa
5.4 3.9 1.7 0.4 setosa 
## 'data.frame':    150 obs. of  5 variables:
##  $ Sepal.Length: num  5.1 4.9 4.7 4.6 5 5.4 4.6 5 4.4 4.9 ...
##  $ Sepal.Width : num  3.5 3 3.2 3.1 3.6 3.9 3.4 3.4 2.9 3.1 ...
##  $ Petal.Length: num  1.4 1.4 1.3 1.5 1.4 1.7 1.4 1.5 1.4 1.5 ...
##  $ Petal.Width : num  0.2 0.2 0.2 0.2 0.2 0.4 0.3 0.2 0.2 0.1 ...
##  $ Species     : Factor w/ 3 levels "setosa","versicolor",..: 1 1 1 1 1 1 1 1 1 1 ...
##  [1] setosa     setosa     setosa     setosa     setosa     setosa    
##  [7] setosa     setosa     setosa     setosa     setosa     setosa    
## [13] setosa     setosa     setosa     setosa     setosa     setosa    
## [19] setosa     setosa     setosa     setosa     setosa     setosa    
## [25] setosa     setosa     setosa     setosa     setosa     setosa    
## [31] setosa     setosa     setosa     setosa     setosa     setosa    
## [37] setosa     setosa     setosa     setosa     setosa     setosa    
## [43] setosa     setosa     setosa     setosa     setosa     setosa    
## [49] setosa     setosa     versicolor versicolor versicolor versicolor
## [55] versicolor versicolor versicolor versicolor versicolor versicolor
## [61] versicolor versicolor versicolor versicolor versicolor versicolor
## [67] versicolor versicolor versicolor versicolor versicolor versicolor
## [73] versicolor versicolor versicolor
##  [ reached getOption("max.print") -- omitted 75 entries ]
## Levels: setosa versicolor virginica

Και ας μετονομάσουμε τα ονόματα των ειδών2:

1.8 Λιγότερα δεκαδικά

Με την ακόλουθη εντολή, θα εμφανίζονται πλέον λιγότερα δεκαδικά:

1.9 Στατιστικά στοιχεία για κάθε taxon

Θέλουμε να μάθουμε την ελάχιστη, την μέγιστη και την μέση τιμή, καθώς και την τυπική απόκλιση3:

## $`Iris setosa`
##   Sepal.Length  Sepal.Width   Petal.Length   Petal.Width              Species  
##  Min.   :4.3   Min.   :2.3   Min.   :1.00   Min.   :0.10   Iris setosa    :50  
##  1st Qu.:4.8   1st Qu.:3.2   1st Qu.:1.40   1st Qu.:0.20   Iris versicolor: 0  
##  Median :5.0   Median :3.4   Median :1.50   Median :0.20   Iris virginica : 0  
##  Mean   :5.0   Mean   :3.4   Mean   :1.46   Mean   :0.25                       
##  3rd Qu.:5.2   3rd Qu.:3.7   3rd Qu.:1.57   3rd Qu.:0.30                       
##  Max.   :5.8   Max.   :4.4   Max.   :1.90   Max.   :0.60                       
## 
## $`Iris versicolor`
##   Sepal.Length  Sepal.Width   Petal.Length  Petal.Width              Species  
##  Min.   :4.9   Min.   :2.0   Min.   :3.0   Min.   :1.00   Iris setosa    : 0  
##  1st Qu.:5.6   1st Qu.:2.5   1st Qu.:4.0   1st Qu.:1.20   Iris versicolor:50  
##  Median :5.9   Median :2.8   Median :4.3   Median :1.30   Iris virginica : 0  
##  Mean   :5.9   Mean   :2.8   Mean   :4.3   Mean   :1.33                       
##  3rd Qu.:6.3   3rd Qu.:3.0   3rd Qu.:4.6   3rd Qu.:1.50                       
##  Max.   :7.0   Max.   :3.4   Max.   :5.1   Max.   :1.80                       
## 
## $`Iris virginica`
##   Sepal.Length  Sepal.Width   Petal.Length  Petal.Width              Species  
##  Min.   :4.9   Min.   :2.2   Min.   :4.5   Min.   :1.40   Iris setosa    : 0  
##  1st Qu.:6.2   1st Qu.:2.8   1st Qu.:5.1   1st Qu.:1.80   Iris versicolor: 0  
##  Median :6.5   Median :3.0   Median :5.6   Median :2.00   Iris virginica :50  
##  Mean   :6.6   Mean   :3.0   Mean   :5.6   Mean   :2.03                       
##  3rd Qu.:6.9   3rd Qu.:3.2   3rd Qu.:5.9   3rd Qu.:2.30                       
##  Max.   :7.9   Max.   :3.8   Max.   :6.9   Max.   :2.50
##              Iris setosa Iris versicolor Iris virginica
## Sepal.Length        0.35            0.52           0.64
## Sepal.Width         0.38            0.31           0.32
## Petal.Length        0.17            0.47           0.55
## Petal.Width         0.11            0.20           0.27

Μπορείτε επίσης να χρησιμοποιήσετε το σύνολο εντολών από την πρώτη εργαστηριακή άσκηση, προκειμένου να δημιουργήσετε τα στατιστικά αυτά δεδομένα. Δοκιμάστε το για το IQR (Interquantile Range).

1.10 Γραφική απεικόνιση των αποτελεσμάτων

Ας αναπαριστήσουμε γραφικά τα δεδομένα μας για κάθε taxon και μια μεταβλητή4:

## [1] "G:/Academic/Lessons/EMB/Labs/Labs/TRUE RMD Files/EKPA/2021"

Κάντε το ίδιο με τις υπόλοιπες μεταβλητές και αλλάξτε το όνομα του y άξονα.

ΠΡΟΣΟΧΗ
Για την εργασία σας, θα χρειαστεί να ελέγξετε εάν οι μεταβλητές σας έχουν κανονική κατανομή και εάν εμφανίζουν μεγάλο συντελεστή συσχέτισης. Ανατρέξτε σε προηγούμενες εργαστηριακές ασκήσεις για να δείτε πώς μπορείτε να το κάνετε αυτό.

2 Κυρίως ανάλυση

Είμαστε πλέον σε θέση να προχωρήσουμε στην κυρίως ανάλυση των μορφομετρικών μας δεδομένων. Είναι το πιο δύσκολο και κρίσιμο βήμα και θα χρειαστεί να είστε εξαιρετικά προσεκτικοί όταν θα τρέξετε την ακόλουθη εντολή για την εργασία σας:

## [1] 0.98

Με την ακόλουθη εντολή βλέπουμε τα αποτελέσματα της ανάλυσης μας:

Lec.lda


## Let's say what we did
##--------------------------------------------------------------------------------##
# Call: lda(Species ~ Sepal.Length + Sepal.Width + Petal.Width + Petal.Length,
# data = iris)
##--------------------------------------------------------------------------------##
## How much of the data lie in each species

## Prior probabilities of groups:

## setosa versicolor virginica

## 0.33 0.33 0.33
##--------------------------------------------------------------------------------##

## Group means:

## Sepal.Length Sepal.Width Petal.Width Petal.Length

## setosa 5.0 3.4 0.25 1.5

## versicolor 5.9 2.8 1.33 4.3

## virginica 6.6 3.0 2.03 5.6
##--------------------------------------------------------------------------------##

## Coefficients of linear discriminants:

## LD1 LD2

## Sepal.Length 0.83 0.024

## Sepal.Width 1.53 2.165

## Petal.Width -2.81 2.839

## Petal.Length -2.20 -0.932

##--------------------------------------------------------------------------------##

## Proportion of trace:

## The most important part

## It shows how many discriminant axes are there and which is the most important
## in order to distinguish the taxa

## The axes are like this: a*x1 + b*x2 We get a & b from the Coefficients of
## linear discriminants


## (here: 0.83*Sepal.Length + 1.53*Sepal.Width - 2.81*Petal.Width -
## 2.20*Petal.Length for the first axis)


## We can see them with: Lec.lda$scaling
##--------------------------------------------------------------------------------##

## LD1 LD2

## 0.9912 0.0088

## LD1 is the most important as it explains 99,12%
##--------------------------------------------------------------------------------##

2.1 Σημαντικότητα του κάθε διαγνωστικού άξονα σε σχέση με τους υπόλοιπους

Θέλουμε η αναλογία των τιμών αυτών, να είναι όσο το δυνατόν μεγαλύτερη, καθώς υποδηλώνουν τη σχετική σημαντικότητα κάθε άξονα. Όσο μεγαλύτερη η τιμή, τόσο πιο διαγνωστικός ο κάθε άξονας:

## [1] 48.6  4.6

2.2 Συσχέτιση των μορφολογικών γνωρισμάτων με τον κάθε διαγνωστικό άξονα

Για την εργασία σας, εάν έχετε περισσότερους από 2 διαγνωστικούς άξονες, θα χρειαστεί να τους προσθέσετε:

##               [,1]
## Sepal.Length -0.79
## Sepal.Width   0.53
## Petal.Length -0.98
## Petal.Width  -0.97
##               [,1]
## Sepal.Length 0.218
## Sepal.Width  0.758
## Petal.Length 0.046
## Petal.Width  0.223

2.3 Πίνακας Ορθής Ταξινόμησης

Σε αυτό το σημείο, θέλουμε να δούμε πόσο καλά διακρίνουν οι μορφολογικοί χαρακτήρες οι οποίοι αναδείχθηκαν ως σημαντικοί από την ανάλυση μας, τα διάφορα taxa μεταξύ τους:

2.4 Σημαντικότητα του κάθε μορφολογικού χαρακτήρα

Τώρα μπορούμε να διαπιστώσουμε πόσο σημαντικός είναι ο εκάστοτε μορφολογικός χαρακτήρας:

## Formula containing included variables: 
## 
## Species ~ Petal.Length + Sepal.Width + Petal.Width + Sepal.Length
## <environment: 0x0000000015026cd0>
## 
## 
## Values calculated in each step of the selection procedure: 
## 
##           vars Wilks.lambda F.statistics.overall p.value.overall
## 1 Petal.Length        0.059                 1180         2.9e-91
## 2  Sepal.Width        0.037                  307        2.8e-103
## 3  Petal.Width        0.025                  258        5.0e-113
## 4 Sepal.Length        0.023                  199        1.4e-112
##   F.statistics.diff p.value.diff
## 1            1180.2      2.9e-91
## 2              43.0      2.0e-15
## 3              34.6      5.1e-13
## 4               4.7      1.0e-02

2.5 Γραφική απεικόνιση των αποτελεσμάτων

Taxon LD1 LD2
Iris setosa 8.1 0.3 
## [1] "Taxon" "LD1"   "LD2"
## [1] "G:/Academic/Lessons/EMB/Labs/Labs/TRUE RMD Files/EKPA/2021"

3 Εργασία για το σπίτι

Σας παρακαλώ θερμά να στείλετε την εργασία σας σε μορφή PDF σε αυτό το e-mail και να απαντήσετε στα ακόλουθα ερωτήματα, αφού έχετε περιγράψει το σετ δεδομένων το οποίο έχετε χρησιμοποιήσει:
1. Πόσα είδη έχει;
2. Πόσα δείγματα έχει το κάθε είδος;
3. Πόσους μορφολογικούς χαρακτήρες περιέχει το σετ δεδομένων σας;
4. Με βάση τα αποτελέσματα σας, να φτιάξετε έναν πίνακα ο οποίος να περιέχει τη μέση, την ελάχιστη και την μέγιστη τιμή για κάθε μορφολογικό χαρακτήρα κάθε είδους, μαζί με την τυπική απόκλιση.
5. Ποιες μεταβλητές δεν εμφανίζουν κανονική κατανομή;
6. Ποιες μεταβλητές εμφανίζουν μεγάλο (\(>0.8\)) συντελεστή συσχέτισης;
7. Εάν υπάρχουν κάποια ζεύγη ανεξάρτητων μεταβλητών με σ.σ. \(>0.8\), θα τις χρησιμοποιήσετε όλες;
8. Πόσοι διαγνωστικοί άξονες εμφανίστηκαν στην ανάλυση σας; Ποιο είναι το ποσοστό της ποικιλότητας που εξηγεί ο καθένας εξ αυτών; Ποιος είναι ο σημαντικότερος και πώς μπορείτε να είστε βέβαιοι για την απάντηση σας;
9. Ποιος μορφολογικός χαρακτήρας σχετίζεται περισσότερο με κάθε εξ αυτών των αξόνων;
10. Δημιουργήστε τον πίνακα ορθής ταξινόμησης. Ποιο είναι το ποσοστό ορθής ταξινόμησης για κάθε taxon το οποίο έχετε στην ανάλυση σας; Υπάρχει ποσοστό αλληλοεπικάλυψης μεταξύ κάποιων taxa; Εάν ναι, που θεωρείτε ότι οφείλεται;
11. Ποιος μορφολογικός χαρακτήρας είναι ο σημαντικότερος για την διάκριση των taxa που έχετε στην ανάλυση σας;
12. Δημιουργήστε ένα γράφημα το οποίο να απεικονίζει γραφικά τα αποτελέσματα της ανάλυσης σας. Τι παρατηρείτε; Συμφωνούν τα αποτελέσματα του πίνακα ορθής ταξινόμησης με το εν λόγω γράφημα;

4 Appendix: R code

##===========================================================================##
## Do NOT run this code
##===========================================================================##
## Global options

library(knitr)
library(rmdformats)

## Global options
options(max.print="75")
opts_chunk$set(echo=TRUE,
                 cache=TRUE,
               prompt=FALSE,
               tidy=TRUE,
               # comment=NA,
               message=FALSE,
               warning=FALSE,
               eval = TRUE)
opts_knit$set(width=75)
##===========================================================================##
pre {
  max-height: 400px;
  overflow-y: auto;
}

pre[class] {
  max-height: 400px;
}
## CRAN repository:
install.packages("easypackages", repos = "http://cran.us.r-project.org") 
library(easypackages)
libraries("rio", "car", "MASS", "psych", "leaps", "ggplot2", "klaR",
 "DiscriMiner","dplyr", "car", "MASS", "corrgram", "tidyverse", "plyr")
getwd() 
#setwd('E:/Academic/Lessons/Labs/Labs')
data(iris)
names(iris) 
head(iris) 
str(iris) 
iris$Species
iris$Species <- revalue(iris$Species, c("setosa"="Iris setosa", 
"virginica"="Iris virginica", "versicolor"="Iris versicolor"))
options(digits = 2)
iris %>% split(.$Species) %>% map(~summary(.))

DiscriMiner::groupStds(iris[,1:4], iris[,5])
## Violin plot
ggplot(iris, aes(Species, Petal.Length)) +
  geom_violin(adjust = 4)  + 
  geom_boxplot(width = 0.1, 
             fill = "black", 
             outlier.colour = NA) +
  stat_summary(fun.y = median, 
             geom = "point",
             fill = "white", 
             shape = 21, 
             size = 2.5) +
  xlab("Taxon") +
  ylab("Petal Length") +
  theme(axis.title = element_text(face = "bold"),
        axis.text.x = element_text(face = "italic", 
                                   colour = "black"))

getwd()  

ggsave("Petal length.png",
       width = 20, 
       height = 20,
       units = "cm", 
       dpi = 600)                           
attach(iris)  ## For your assignment, change that to mydata

iris <- iris[apply(iris,1,function(row) all(!is.na(row))),] ## Do the same here


set.seed(2934) ## Change that if you want to whichever number you like

cv.grp <- sample(1:10, nrow(iris), 
                 replace = TRUE) ## Remember to change the name

cv.lda <- c() ##A function that splits our data to 2 independent subsets, one to
              ## train our model and one to test how well our model behaves

for(i in 1:10) {
  train.ind <- which(cv.grp!=i)
  test.ind <- which(cv.grp==i)
  model.train <- lda(Species ~ Sepal.Length +## Change the variable names
                                             ## according to the ones that 
                                             ## appear in your data

  Sepal.Width + Petal.Width + Petal.Length,  ## Do the same here
  data = iris, subset = train.ind)           ## Change to mydata
  cv.lda[test.ind] <- predict(model.train,
  newdata = iris[test.ind,])$class             ## Do the same here
 }

mean(cv.lda==as.numeric(iris$Species))       ## And here

Lec.lda <- lda(Species ~ Sepal.Length +      ## This is the test point
Sepal.Width + Petal.Width + Petal.Length,    ## Change the variable names

data = iris) ## Change to mydata


plot(Lec.lda,col = c("green","blue","red", "black", "orange") 
     
[iris$Species], pch = 19) ## Change to mydata


Lec.lda 


## Let's say what we did
##--------------------------------------------------------------------------------##
# Call: 
# lda(Species ~ Sepal.Length + Sepal.Width + Petal.Width + Petal.Length, 
#    data = iris)
##--------------------------------------------------------------------------------##
## How much of the data lie in each species

## Prior probabilities of groups:  

##    setosa versicolor  virginica 

##     0.33       0.33       0.33 
##--------------------------------------------------------------------------------##

## Group means:
            
           ## Sepal.Length Sepal.Width Petal.Width Petal.Length

## setosa            5.0         3.4        0.25          1.5

## versicolor        5.9         2.8        1.33          4.3

## virginica         6.6         3.0        2.03          5.6
##--------------------------------------------------------------------------------##

## Coefficients of linear discriminants:  

##                LD1    LD2

## Sepal.Length  0.83  0.024

## Sepal.Width   1.53  2.165

##Petal.Width  -2.81  2.839

##Petal.Length -2.20 -0.932

##--------------------------------------------------------------------------------##

## Proportion of trace:  

                       ## The most important part

                       ## It shows how many discriminant axes are there and 
                       ## which is the most important in order to distinguish 
                       ## the taxa 

                       ## The axes are like this: a*x1 + b*x2 
                       ## We get a & b from the Coefficients of linear  
                       ## discriminants

                       ##(here: 0.83*Sepal.Length + 1.53*Sepal.Width 
                       ## - 2.81*Petal.Width - 2.20*Petal.Length for the first axis)

                       ## We can see them with: Lec.lda$scaling
##--------------------------------------------------------------------------------##

##   LD1    LD2 

## 0.9912 0.0088   

## LD1 is the most important as it explains 99,12%
##--------------------------------------------------------------------------------##
SVD <- Lec.lda$svd
SVD
LD1 <- predict(Lec.lda)$x[,1]
LD2 <- predict(Lec.lda)$x[,2]

## We will see which variable is the most important for each axis
## We are interested in the absolute value

cor(iris[,1:4],LD1) ## Change to mydata
                    ## Change the number of columns

cor(iris[,1:4],LD2) ## Do the same
x <- Lec.lda
y <- predict(x, iris) ## Change to mydata

detach(iris)        ## Change to mydata - DON'T FORGET TO DO THAT

round(100*errormatrix(iris$Species, y$class, relative=T), 1) ## Change to mydata

## Change to mydata

## The smaller the Wilk's lambda, the better

klaR::greedy.wilks(Species ~ Sepal.Length + Sepal.Width +
               Petal.Width + Petal.Length, 
             data = iris, 
             niveau = 0.1) 
                                                       

s <- predict(Lec.lda)$x
s <- as.data.frame(s)
Taxon <- factor(iris$Species) ## ## Change to mydata 

s <- cbind(Taxon, s)
head(s,1)
names(s)

getwd() 

ggplot(s, aes(x = LD1)) + 
  geom_histogram(fill = "white",
                 colour = "black")  + 
  facet_grid(Taxon ~ .)                         

ggsave("Iris LD1.png", width = 20, height = 20,
       units = "cm", dpi = 600) 

ggplot(s, aes(x = LD2)) + 
  geom_histogram(fill = "white", 
                 colour = "black")  +
  facet_grid(Taxon ~ .)                         
                                            
ggsave("Iris LD2.png", 
       width = 20, 
       height = 20,
       units = "cm", dpi = 600)                       


a <- s$LD1
b <- s$LD2
TAXON <- iris$Species ## ## Change to mydata

g <- data.frame(TAXON, a, b) ## Create a dataframe that contains all the above

## Another dataframe containing the centroids for each species and each
## discriminant axis
centroids <- aggregate(cbind(a,b) ~ TAXON, 
                       g, 
                       mean) 

ggplot(g,aes(a, b, 
             color = TAXON, 
             shape = TAXON))  +
  geom_point(size = 3)  +
  stat_ellipse(geom = "polygon", 
               alpha = 1/2, 
               aes(fill = TAXON))  +
  xlab("LD1") +
  ylab("LD2") +
  theme(axis.title = element_text(face = "bold"))  +
  stat_summary(data = centroids,
               fun.y = "mean",
               geom = "point",
               shape = 22, 
               size = 4,
               colour = "black", 
               aes(fill = TAXON))  +
  
  ## Shape of each species 
  
  ## If you have more than 3 species you must add some values (range: 10-20)
  scale_shape_manual(values = c(15, 16, 17))  
                                              

ggsave("Iris ellipses.png", 
       width = 20,
       height = 20,
       units = "cm", 
       dpi = 600)                            

##===========================================================================##
## Do NOT run this code
##===========================================================================##
labs = knitr::all_labels(echo = TRUE)
##===========================================================================##

  1. Hint: Ανατρέξτε σε προηγούμενες εργαστηριακές ασκήσεις για να δείτε πώς γίνεται ή διαβάστε το manual της εν λόγω βιβλιοθήκης: readxl↩

  2. Δεν θα χρειαστεί να το κάνετε αυτό εσείς για την εργασία σας↩

  3. Για την εργασία σας, για να βρείτε την τυπική απόκλιση, θα χρειαστεί να τρέξετε την ακόλουθη εντολή πρώτα: mydata.sd <- as.data.frame(mydata) - εννοείται στην περίπτωση που έχετε ονομάσει το αρχείο σας mydata↩

  4. Για να κάνετε αυτό το βήμα για την εργασία σας, θα πρέπει να αλλάξετε το όνομα του αρχείου, καθώς και τα ονόματα των μεταβλητών - Εάν έχετε απορίες, τώρα είναι η στιγμή να τις εκφράσετε↩

 

Delivered to you by Dr. Kostas Kougioumoutzis

kkougiou@aua.gr